Temperature dependence of density profiles for a cloud of non-interacting fermions 
moving inside a harmonic trap in one dimension 
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We extend to finite temperature a Green's function method that was previously proposed to 
evaluate ground-state properties of mesoscopic clouds of non-interacting fermions moving under 
harmonic confinement in one dimension. By calculations of the particle and kinetic energy density 
profiles we illustrate the role of thermal excitations in smoothing out the quantum shell structure of 
the cloud and in spreading the particle spill-out from quantum tunnel at the edges. We also discuss 
^v^j ' the approach of the exact density profiles to the predictions of a semiclassical model often used in 

the theory of confined atomic gases at finite temperature. 
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I. INTRODUCTION 



The achievement of Bose-Einstein condensation in ultracold gases of bosonic atoms has given great impulse to 
the study of dilute quantum gases inside magnetic or optical traps Q . Similar cooling and trapping techniques are 
being used to drive gases of fermionic alkali atoms into the quantum degeneracy regime Q|. It is also possible to 
experimentally generate and study strongly anisotropic atomic fluids effectively approaching dimensionality D = 2 
' O • or D = 1, at very low temperature and with very high purity, inside magnetic traps where the transverse and the 
longitudinal confinement are vastly different ||. 

In the relevant conditions of temperature and dilution the atomic interactions become negligible in a gas of fermions 
placed in a single Zeeman sublevel inside a magnetic trap Q . One can obtain in this way a close laboratory realization 
of an inhomogeneous, non-interacting Fermi gas, which has special significance in regard to the foundations of density 
functional theory and to tests of the Thomas- Fermi approximation ||. In the quasi-onedimensional (ID) case the 
ground state and the excitation spectrum of such an ideal gas of spin-polarized (or " spinless" ) fermions can be mapped 
into those of a gas of hard-core impenetrable bosons ||. The latter is known as the Tonks gas, from the early work of 
Tonks JtJ on the equation of state of hard-object adsorbates. Advances in atom waveguide technology, with potential 
l/^ . applications to atom interferometry and integrated atom optics, especially motivate theoretical studies of dilute gases 
CO ' in a regime where the dynamics becomes essentially ID M. 

The wave functions of a cloud of spinless fermions under harmonic confinement are, of course, Slater determinants 
of harmonic-oscillator single-particle orbitals. Their representation in terms of Hermitc polynomials has, however, 
limited usefulness for numerical calculations on mesoscopic clouds || . Brack and van Zyl jlO) have developed a more 
powerful method for non-interacting fermions occupying a set of closed shells under isotropic harmonic confinement in 
D dimensions, leading to analytical expressions for the particle and kinetic energy densities at temperature T = in 
■ terms of Laguerre polynomials. These expressions are especially useful for numerical calculations on isotropic systems 
in D = 2 and D — 3. A Green's function method, which altogether avoids the use of wave functions in favor of the 
matrix elements of the position and momentum operators, has been developed for similar purposes fTl|| and extended 
CD ■ to calculations of the pair distribution function in the ground state for the ID ideal Fermi gas ]l2j . This method has 
, also been extended to the ground state of ideal Fermi gases under harmonic confinement of arbitrary anisotropy in 
higher dimensionalities [jl3|l . 

The purpose of the present paper is to extend the Green's function method to an inhomogeneous ID cloud of non- 
interacting fermions at finite temperature and to illustrate its usefulness by numerical calculations of the particle and 
kinetic energy density profiles as functions of temperature in the case of harmonic confinement. Analytical expressions 
have been derived by Wang |b|| for the same system at high temperature, where the chemical potential is lower than 
the single-particle ground-state energy. The case of non-interacting fermions under 3D harmonic confinement has 
been treated by Schneider and Wallis |fj| through the use of Laguerre polynomials. The emphasis of our numerical 
calculations will be to illustrate how the characteristic quantum features of the fermion cloud, i.e. its shell structure 
and the spill-out of particles at the boundaries beyond the Thomas-Fermi radius, evolve with increasing temperature 
as a semiclassical regime is being approached. 
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The paper is organized as follows. Section II reports some essential definitions, starting from the one-body density 
matrix at T ^ 0, and gives the expressions of the density profiles in the semiclassical regime. Sections III and IV 
present the essential details of our method and our numerical results, respectively. A brief summary concludes the 
paper in Section V. Although we refer through the paper to the results as being appropriate to spinless fermions in 
ID, they are equally valid for the Tonks gas of impenetrable bosons. 

II. ESSENTIAL DEFINITIONS AND THE SEMICLASSICAL LIMIT 

The generalized grand-canonical density matrix for fermions can be written as 



D(xi,x;/3,n) = ^ — 



1 



exp [j3{Ei - fi)] 



■ip(x — x\ ) 



ipi (xi). 



(1) 



Here (3 = l/fc^T and /i is the chemical potential, while ipi and Ei are the single-particle orbitals and the corresponding 
energy eigenvalues. The zero temperature limit of Eq. (]]]) leads to the Dirac density matrix for N ideal spinless 
fermions. 



(2) 



The particle density profile n(x) of the gas at temperature T and chemical potential /z is the zero-order moment of 
the matrix D(x,x\), 



n(x) = D(xi,x;P,fi)\ Xl= 



:(lpi \S(x - x)\ Ipi). 



,_, ■ ■ exp [/3{Ei - (j,)} 
The kinetic pressure P(x) is given by a specific second-order moment of D(x,xi) jn 

m dx 2 

1 

2m 



(3) 



P(x) = - — -^D(xi,x;P,^)\ Xl=x 



^ OO 

— Y- 

i=i 



- exp [(3{Ei - /x)] 



(ipi \p 2 5(x — x) + S(x — x)p 2 \ ipi). 



(4) 



The kinetic pressure P(x) is twice the kinetic energy density of the fermion cloud. In these equations p and x are the 
momentum and position operators, and the chemical potential zx is determined from the condition 



n(x; /3, fxjdx = N 



(5) 



where N is the average number of particles. 

In the semiclassical regime the particle density and the kinetic pressure of an ideal Fermi gas confined by a ID 
potential V(x) can be calculated by treating the energy levels as a continuum. The appropriate condition of validity 
is that the level spacing AE is sufficiently smaller than the thermal energy kgT. The appropriate expressions in the 
grand-canonical ensemble are 



n sc (x) 



dE S^h 5 { E -h- V{X) ) l + exp^-M)] 
dp 



2nh 



exp 



0[^ + V(x)-, 



1 



and 



dp p 2 



— <i exp 



P3[^ + V(x)-» 



(6) 



(7) 



The chemical potential is again determined by normalization to the average number of fermions. 
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III. THE GREEN'S FUNCTION METHOD 



Equations (|^) and (0) can be rewritten in terms of the Green's function G(x) = \im £ _, + (x — x + ie) 1 in coordinate 



space, 



1 00 1 

n ( x ) = lim Im — — — ^{if>i\d(x)\i/)i 

7r e-+o+ f— ' 1+exp \B[Ei — u)\ 



_ { ,exp[/3(£;- M )] 
— lim ImTr (T ■ G(x)) (8) 

7T e^0+ 



and 



1 00 1 "2 

P(x) = lim Im V — -{^ \—G(x)\ ipi) 

v; 7re^o+ ^ l + exp[/3(^ -/*)] "m v 71 ^' 

= — lim ImTr (T ■ f ■ Gix)). (9) 

7rm e^o+ 



Here we have used a matrix formalism, by introducing the temperature matrix T whose diagonal elements [T]^ = 
1/{1 + exp \P(Ei — /i)]}^ 1 contain the statistical Fermi factors, while the off-diagonal elements \T\i,j with i ^ j are 
null. 

The evaluation of the expressions in Eqs. (||) and (g) is carried out by an immediate extension of the procedure 
developed in pd| ] for the case of zero temperature. The trace of a generic matrix A is connected to the elements of 
the inverse matrix A^ 1 by the relation 



TrA = ^ [lndet^+AI)])^. (10) 



We get from Eqs. (§) and @ 



Id 

n(x) = lim Im— [lndetfx — x + XT + ie)] ,_ n (11) 

7T e^0+ 5 A /JiA-u 



and 



1 d 

P(x) = lim Im— [In detfa; - x + ie + XT ■ p 2 )} , n . (12) 

nm e^o+ oX IA_U 

In the specific case of harmonic confinement we make use of the representation of the position and the momentum 
operators in the basis of the eigenstates of the harmonic oscillator with energy E n — (n — l/2)hw: that is, x = 
(a + a^)/y/2 and p = i(aJ — a)/ y/2 with a \ ip n ) = y/n — 1 1 ip n -\) and \ if> n ) = y/n \ ip n ^i). The tridiagonal form 
of the matrices representing x and p allows us to express the determinants in Eqs. ( pi] ) and (|l^) by a recursive 
algorithm [pd|Jl7| as products of an infinite number of matrices having dimension lxl for the particle density and 
2 x 2 for the kinetic pressure. If we write the matrix JC V = x — XT ■ p v in the tridiagonal form 



K." 



(A x £1,2 

62,1 A2 B 2 ,3 J (13) 



where the indices refer to blocs of dimension 1 x 1 if v = and 2 x 2 if v = 2, we obtain 

oo 

det(x - K! J + ie) = TJ det(sc - Aj + ie). (14) 
i=i 

Here Ai = A\ while the renormalized blocs Aj for j > 1 are given by the recursive formula 

Aj = Aj + Bjj-x(x - Aj-i + iey-Bj-ij. (15) 

As in the calculation of density profiles at T = jll| the evaluation of the determinant in Eq. ( |i"4| ) converges quite 
rapidly, yielding accurate results for mesoscopic fermion clouds with a moderate amount of numerical effort. 
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IV. NUMERICAL RESULTS 



We report in this section some numerical results for the particle density and the kinetic pressure that we have 
obtained by the Green's function method. Our main purpose is to illustrate the approach of the profiles to their 
semiclassical expressions and we consider for clarity clouds containing rather small average numbers of particles, i.e. 
N = 4 and N = 20. Some sharpening of the structures in the profiles should be expected in a canonical ensemble 
viewpoint, where fluctuations in the particle number are suppressed. 

In the calculations reported below we have approximated the determinant in Eq. (|lj) by the product of its first 
M terms with M — 10 7 as an upper bound. This choice allows for up to 10 7 thermally excited states and the 
corresponding value of the spectral resolution parameter is chosen as e = 10~ 3 . These two choices evidently limit to 
some minor extent the accuracy of our numerical results relatively to the formally exact theory given in Sect. III. 
Of course, for such small numbers of particles one may also approach the same calculations by using directly the 
expression of the harmonic-oscillator orbitals in Eq. (0). However, the routines that are available to evaluate Hermite 
polynomials also involve numerical approximations which are increasingly severe for polynomials of high degree, as 
needed to describe the gas at finite temperature. 

The first step in our calculations is to evaluate the chemical potential at fixed N by means of a self-consistent 
solution of Eqs. (||) and (||) |Q. The resulting semiclassical value of \x is used in Eqs. ( |ll"| ) and (12) for the evaluation 
of the particle density and kinetic pressure profiles. The procedure is justified a posteriori by the fact that integration 
of the spatial density n(x) obtained at this stage reproduces the correct average number N of fermions. 

In Figures 1 and 2 we report the particle density profiles at various temperatures for N — 4 and N — 20, respectively. 
The corresponding profiles of kinetic pressure are given in Figures 3 and 4. In each Figure the left-hand panel shows 
how the quantum effects, which consist of oscillations from shell structure in the profiles and in particle spill-out 
and negative kinetic pressure at the edges, are washed away by thermal excitations. The right-hand panels show 
instead how the profiles approach the semiclassical regime: the latter already holds to a good approximation when 
ksT ~ hui, although some trace of a negative kinetic pressure in the spill-out region remains in the "exact" profile at 
this temperature. 



V. CONCLUSIONS 



In conclusion, in this paper we have extended to finite temperature the Green's function method that was previously 
developed for the formally exact evaluation of ground-state properties of an inhomogeneous ideal Fermi gas in ID. 
The method dispenses with knowledge of the single-particle orbitals in favour of the single-particle energy levels and of 
the matrix elements of the position and momentum operators. Our numerical applications to clouds containing small 
numbers of particles under harmonic confinement have shown that the commonly used semiclassical approach to a 
thermally excited fermion cloud yields (i) accurate values of the chemical potential except at the lowest temperatures, 
and (ii) a reliable account of particle and kinetic energy density profiles once their main quantum features become 
essentialy washed out at UbT ~ fku. 

Two further comments are in order. Firstly, the profiles shown in Figures 1 and 2 also describe (aside from an 
immediate rescaling of the units used on the axes) the momentum distribution of the fermion cloud under harmonic 
confinement. Secondly, the same profiles also refer to the Tonks gas of impenetrable bosons in ID, as was explicitly 
shown by Yang and Yang p9| ]. 
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FIG. 1. Particle density profile for TV = 4 harmonically confined fermions at various values of the temperature. Left panel: 
"exact" profiles at T = (solid curve, calculated by the method of Ref. jnj) and at T = 0.2hu/kB (dashed curve); the other 
curves refer to ksT/hui = 0.5, 1.0, 2.0 and 3.0, in order of decreasing peak height. Positions are in units of the harmonic 
oscillator length a^o = \J h/(mcu) and the particle density is in units of . The right-hand panel reports again the profiles at 
kBT/hui = 0.2, 1.0 and 3.0, together with those calculated in the semiclassical approximation (light lines); for the latter two 
values of fcsT/fitj the "exact" profiles are hardly distinguishable from the semiclassical ones. 
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FIG. 2. The same as in Figure nl, for iV = 20 harmonically confined fermions. 
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FIG. 3. Kinetic pressure profile for N = 4 harmonically confined fermions at the same values of the temperature as in 
Figure [I]. Left panel: "exact" profiles; right panel: comparison with the semiclassical profiles for hsT/hcd — 0.2, 1.0 and 3.0 
(light lines). Positions are in units of a^o and the kinetic pressure is in units of ftwa^. 
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